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Abstract 

A complete numerical implementation, in both singlet and non-singlet sectors, 
of a very elegant method to solve the QCD Evolution equations, due to Furmanski 
and Petronzio, is presented. The algorithm is directly implemented in x-space by a 
Laguerre expansion of the parton distributions. All the leading-twist distributions 
are evolved: longitudinally polarized, transversely polarized and unpolarized, to 
NLO accuracy. The expansion is optimal at finite x, up to reasonably small x- 
values (x ~ 10~ 3 ), below which the convergence of the expansion slows down. The 
polarized evolution is smoother, due to the less singular structure of the polarized 
DGLAP kernels at small-x. In the region of fast convergence, which covers most 
of the usual perturbative applications, high numerical accuracy is achieved by 
expanding over a set of approximately 30 polynomials, with a very modest running 
time. 
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Program Summary 

Title of the Programs: nsunpol, snunpol, nslong, snlong, trans 

Computer: 
Sun 19 

Operating system: Unix 

Programming language used: FORTRAN 77 

Peripherals used: Laser Printer 

Number of lines in distributed program: 4260 

Keywords: Structure function, polarized parton distribution, Q 2 evolution, Laguerre ex- 
pansions, numerical solution. 

Nature of physical problem: 

The programs provided here solve the DGLAP evolution equations, with next-to-leading 
order a s effects taken to account, for unpolarized, longitudinally polarized and trans- 
versely polarized parton distributions. 

Method of solution: 

The method developed by Furmanski and Petronzio is used. The kernel P(x) of the 
DGLAP integrodifferential equations and the evolution operators E(t,x) are expanded 
in Laguerre polynomials. 

Typical running time: 

About 5 seconds for the transverse polarization case and 30 minutes for the longitudinal 
polarization and for the unpolarized. 

LONG WRITE-UP 
1 Introduction 

The study of the spin structure of the proton is a fascinating aspect of the theory of the 
strong interactions. Parton distributions in the nucleon tell us about the structure of 
fundamental observables such as spin, parton densities and correlations among partons, 
in a light-cone framework. In a second quantized theory they emerge as matrix elements 
of nonlocal operators at light-like separations. At low energy they can be described by 
valence quark models and the impact of the QCD evolution (or Renormalization Group 
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Evolution, RGE) on the shape of these distributions can be performed within the parton 
model. 

It has become more and more common in the high energy literature to describe the 
evolution of parton distributions by starting from a low energy input. 

There is widespread interest in the analysis of the effect of the QCD evolution es- 
pecially for those leading twist distributions (such as the transverse spin distribution) 
which have not yet been measured. For instance, in ref. ||, the authors analize the 
impact of the evolution on transverse spin distributions derived at low Q from the Isgur- 
Karl quark model 0, and these predictions can be used direclty -at very high energy- 
for other predictions, such as in the study of the transverse spin dependence of the Drell 
Yan process. A similar analysis can be done for the chiral quark model ||. 

On the numerical side, various methods have been presented in the literature [FTO 



which all try to solve the DGLAP equations by iterating the evolution over infinitesimal 
steps in the fractional momentum x. Another technique is based on the use of Mellin 
moments and on their inversion. In general, moments are equivalent to finite -rather than 
infinitesimal- discretizations of the integro-differential equations. In the case of the Mellin 
moments the inversion (to x-space) is the real difficult and time consuming part of the 
method. 

A different, and very elegant implementation of methods based on finite discretiza- 
tions of RGE's for QCD was formulated long ago by Furmanski and Petronzio 0. Their 
method uses the Laguerre expansion of the initial distributions and of the kernel of the 
evolution equations arrested to an arbitrary order n. The algorithm that they provide 
defines the structure of the moments recursively in terms of some initial conditions. In 
the non-singlet sector there have been attempts to apply the method both to leading and 



to next to leading order |T8|, [Tl| , and to leading order in the singlet sector as well [|T| 



However, a complete implementation and numerical documentation of this algorithm is 
still missing. We also mention that our interest in this method has aroused from our 
search for the numerical implementations of more involved evolutions equations, such as 
those describing the dynamics of Compton scattering in the deeply virtual limit (DVCS) 
||. In this latter case, the RGE evolution is a continuous interpolation between 2 lim- 
its: the DGLAP (or Altarelli Parisi) evolution and the Efremov-Radyushkin-Brodsky- 
Lepage evolution (ERBL) M. We believe that a complete numerical understanding of 
these "non-diagonal" RGE's and their robust numerical implementation requires finite 
step integrations. We hope to get back to this issue in the near future, and for the rest 
of this paper just focus on the usual DGLAP evolution. Here we have implemented 
the evolution of all the leading twist parton distributions, using the kernels calculated 
by various authors [fj, |3|, [13], A complete list of these kernels can be found in the 
Appendix. 
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2 The Laguerre expansion 



The Laguerre method for the numerical solution of the evolution equations is due to 
Furmanski and Petronzio. In this section we briefly outline the method, which is fast 
converging at intermediate values of x and any Q 2 value. At very small-x (approximately 
10" 3 ), the Laguerre expansion suffers from numerical instabilities, due to the growth of 
the moments. We start by defining our notation and other conventions. 
The two-loop running of the coupling constant is defined by 

a( Ql)_2 1 fl ftlnMcy/A-) 1 \ 



where 



2tt /3 ln(Q 2 /A 2 ) V Po hi(Q 2 /A 2 ) Tn 2 (Q 2 /A 2 )' 



Pi = y°g - Y C ° nf ~ 2Cpnf (2) 



where 



N 2 - 1 1 
C G = N, C F = -^-, T R = - (3) 

and N is the number of colours. 

The solution for the running coupling is given by 

with a(Qo) = a(0), and Qo denoting the initial scale at which the evolution starts. The 
evolution equations are of the form 



Q 2 Jj-M x >Q 2 ) = ^-P(-)(x,a(Q 2 ))^ X i(x,Q 2 ) 



(5) 



with 



Xi (x, Q 2 ) = q^\x, Q 2 ) - —q^(x, Q 2 ) (6) 
for the non-singlet distributions and 



n 2 _±_ ( q (+ \x,Q 2 ) \ _ ( P qq (x,Q 2 ) P qg (x,Q 2 ) \ ( qM(x,Q 2 ] 
W dQ 2 { G(x,Q 2 ) j {P gq (x,Q 2 ) P gg (x,Q 2 ) r{ G(x,Q 2 ) 



(7) 
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for the singlet sector. 

We have defined, as usual 



n f 

q\~ ] = a - ft, q\ +) = ft + ft, <? (+) = s = J2 <?i +) 

»=i 

We introduce the evolution variable 



which replaces Q 2 . The evolution equations are then rewritten in the form 



^qt\t,x) = (pW( x ) + ^R H ( x ) + ..)®q(-\t,x) (10) 



dt 



Q 2 ^Xi(x,Q 2 ) = [P^(x) + ^R {+) (x))®Ux,Q 2 ), (ID 



(it' 

_d / gW(t,i) 
I G(x,t) 



P(°)(x) + ^i? (+) (x)^®^(x,g 2 



In the new variable t, the kernels of the evolution take the form 



= P (1, (x)-|p (0) (x). (13) 



Equations (|10j) and ([□]) are solved independently for the variables q\ ^ and Xi respec- 
tively. Finally, the solution of eq. ( |12|) (or the singlet equation) is substitued into 
Xi in order to obtain q\ + ^ ■ The equations can be written down in terms of two singlet 
evolution operators E±(t,x) and initial conditions q±(x,t = 0) = g±(x) as 

j t E ± = P ± ®E ± , (14) 

whose solutions are given by 

qt\t,x) = £(_) ® gH 

= £(+) ®Xi(aO- ( 15 ) 
The singlet evolution for the matrix operator E(t,x) 



EpF Epc 

Egf Eqg 



(16) 
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is solved similarly as 



d 4 = P®E 
at 



(17) 



(19) 

The method of Furmanski and Petronzio requires an expansion of the splitting func- 
tions and of the parton distributions in the basis of the Laguerre polynomials 



MV) = £^](-1>^ PO) 
which satisfies the property of closure under a convolution 

L n (y) <g> L m {y) = L n+m (y) - L n+m+1 (y). (21) 

In order to improve the small-x behaviour of the algorithm, from now on, the evolution 
is applied to the modified kernel xP(x), which, for simplicity, is still denoted as in all the 
equations above, i.e. by P(x). At a second step, the < x < 1 interval is mapped into 
an infinite interval < y < oo by a change of variable x = e~ y and all the integrations 
are performed in this last interval. We start from the non-singlet case by defining the 
Laguerre expansion of the kernels and the corresponding (Laguerre) moments to lowest 
order 

oo 

P$\V) = E^n(l/), 
n=0 

/'DC 

P(°) = / dye-VLM P(°)(y) (22) 
Jo 



and to NLO 



R(y) = J^R n L n (y). (23) 



n=0 

One defines also the difference of moments 



p (0 = ^0)_^(0) (p (o, =0) 
r, = Ri-Ri-! R-! = 0. (24) 
A similar expansion is set up for the evolution operators E(t, y) 
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E^(t,y) = J2E^(t)L n (y) 

n=0 

oo 

E{t,y) = J2 E n(t)Ln{y), (25) 



n=0 



where all the information on the t evolution is contained in the moments E n (t). The 
solution to NLO is expressed as M 



Mt) = E<°Ht) - ±mpm E m {t) , (26) 



2 a(t)-a(0) (1) 

where 

E^(t) = e p o'^^- (27) 

fc=0 

n 

E^(t) = Y,r n -^°\t), (28) 

(29) 

and the coefficients A^> are determined recursively from the moments of the lowest order 
kernel P(°> 

4 0) = 1 

4 fc+1) = EAf (* = 0,l,2,..,n-l). (30) 

i=k 

In the singlet case one proceeds in a similar way. The solution is expressed in terms 
of a 2-by-2 matrix operator 

oo 

E^(t,y) = ^E^(t)L n (y). (31) 

n=0 

The solution (at leading order) is written down in terms of 2 projection matrices and 
one eigenvalue (A) of the (matrix) kernel 



61 A 



\P (0 \ e 2 = I + Al) , (32) 

(33) 
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where 

A = -{\c F + 2 -n f T R ), (34) 

in the form 

The recursion relations which allow to build A^> and B^ are solved in two steps as 
follows. One solves first for two sets of matrices and by the relations 



4 0) = o 

n-l 



(0) (k) 



i=k 

bf = 

n—1 



= -Ae# + E^A W , (36) 

i=l 

which are used to construct the matrices and 



4 0) = e 2 - 1 (e lfl W - (-l)"e#) 

5i 0) = ei + ^ (e^) - (-ire#) ■ (37) 



These matrices are then input in the recursion relations 



4 0) = e 2 B<°> = ei 

n-l 

JO) n (fc) 



i=k 



(38) 



= -Ae^ + E^W (39) 



witn n > and k — 0, 1, — 1, which generates the coefficients of the matrix-valued 
operator (i.e. the leading order solution). The NLO part of the evolution is obtained 
from 
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EW(t,y) = J2E£\t)L n (y), (40) 

n=0 

with 

E£\t) = E£\t)-2E^ 1 (t) + E^ 2 (t) (41) 

where 

E0)(t) = f dre'^ 2 E (0) (t - ^RjE^ (r)5(n -i-j-k). (42) 

The expressions of E^ and E^ are inserted into eq. (p6j) thereby providing a complete 
NLO solution of the singlet sector. 



3 The polarized and the unpolarized evolution 

The implementation of the polarized and of the unpolarized evolution is performed in the 
MS scheme, which is by now standard in most of the high energy physics applications. 
In the unpolarized case, we introduce valence quark distributions qv(%,Qo) and gluon 
distributions G(x, Qq) at the input scale Qo, taken from the CTEQ parametrization [^5] 



q(x) = A x Al (l -x) M {l + A 3 x Ai ). (43) 

Specifically 



xu v (x) = 1.344a; - 501 (1 - x) 3 ' 689 [l + 6.042a; - 873 ] 
xd v (x) = 0.640a; - 501 (1 - a;) 4 - 247 [l + 2.690a; - 333 ] 
xG(x) = 1.123aT a206 (l -x) 4 - 673 [l + 4.269a; 1 - 508 ] (44) 

and an asymmetric sea contribution 

X q(x) = i[0.255x-°- 143 (l - x) 8 - 041 (l + 6.112a;) T 0.071x°- 501 (l - x) 8 041 ], (45) 

where the (— ) holds for the u and the (+) for the d flavors. The set accounts for a 
u, d flavour asymmetry, and the sea quark contribution is parameterized by 



xs(x) = [0.064a;- - 143 (l-a;) 8U41 (l + 6.112a;)]. (46) 



In the polarized case we have chosen the first set of ref. |£0| which is of the functional 
form ABx c (l - x) D (l + Ex + Fy/x) 
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xAu v = 0.918* 1.365 *a; a512 (l-x) 3 - 96 (l + 11.65x-4.6v / ^) 
xAd v = -0.339 * 3.849x a78 (l-x) 4 - 96 (l + 7.81x-3.48v^) 
xAG = 1.71 * 3.099x a724 (l-x) 5n (l + 0.0x + 0.0^). (47) 

For the (flavour symmetric) sea contribution we have set Am = Ad = As with 

As = -0.06 * 18.521x a724 (l - x) 14 ' 4 (l + 4.63x - A.96y/x). (48) 



Notice that the evolution in the MS requires some care, depending upon the way the 
subtraction of the collinear singularities in the coefficient functions (hard scatterings) 
is performed. The (non singlet) hard scatterings, in fact, do not conserve helicities in 
the annihilation channels. In a "traditional" MS scheme, one has to keep both these 
helicity violating terms of the coefficient functions and has add to the polarized non- 
singlet kernels of the Appendix some additional terms proportional to (1 — x) to NLO [[15 



IT , [i~3|| As a result of this, both the singlet and the non-singlet evolution are affected to 
NLO. However, one can factorize out of the coefficient functions these spuriuos (helicity- 
violating) terms and absorb them into the renormalization of the parton distributions. 
As a result of this procedure, the NLO kernels of the evolution turn out to be exactly 
those defined in the Appendix and the hard-scatterings are helicity preserving. 



4 The evolution of the transverse spin distribution 



The first identification of the transverse spin distribution is due to Ralston and Soper J23] 
in their study of the factorization formula for the Drell Yan cross section. The parton 



interpretation of this distribution has been discussed in various papers [Tj^ [3] , in which 
its behaviour as a leading twist distribution has been pointed out. 

It appears in the double transverse spin asymmetry Att for the Drell Yan process 



23, H 



A TT = sin2 cos 2( P e 2 A T gi(xi) A T qi(x 2 ) - 4g . 
TT l + cos 2 6> Eie 2 gi(a:i)gi(x2) 

In eq. fl49|) the angles 9 and <fi are the polar and the azimuthal angles of the momentum 
of one of the two leptons, measured with respect to the beam (6) and to the photon 
polarization directions ((f)). The asymmetry disappears if the momenta of the two leptons 
are both integrated over. In the parton model, A^(g) is interpreted as the probability of 
finding a quark with spin polarized along the transverse spin of a transversely polarized 
proton minus the probability to find it polarized oppositely. This distribution does not 
couple to gluons and is therefore, purely non-singlet. The LO anomalous dimensions 



of this distribution have been calculated in fl6 |, while the NLO corrections have been 
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derived by various authors using both the Operator Product Expansion and the 
method of ref. [12|, extended to the polarized case. 

Since the gluons don't couple in the evolution, the equation is simply written as 







-A T q±{x,Q' 1 



a s (Q 2 



-A T P q ±(x)®A T q±(x,Q 2 ). 



d\nQ 2 ' 2tt 

As in the previous sections, we have set Arq± = A T q ± Arq. 



(50) 



5 Description of the program 

The first step in the calculation is the computation of Laguerre moments -or expansion 
coefficients- of the kernel which are given by an explicit recursion relation. 

The solutions of the integral equations are obtained by first discretizing the integrals 
in the form 

n 

dqf(q) >^Wif(qi), (51) 

i=i 

where Wi are integration weights for the grid-points qi. Various sets of Gauss-Legendre 
grid-points are provided (by the GAUSS and the LEGENDRE subroutines) in the inter- 
val (—1,1). In order to map the grid points and the weights from the interval (—1, 1) 
to the interval (0, oo), one can use various mappings. The possible types of mapping 
provided in the code are 

(i) MAPNO=l: 

V(x) = R mm + ,_ T , 2 { lt X) _ R y (52) 

1 dj -|- z,/ \lt>max ^min) 

(ii) MAPNO=2: 

y{x) — Rmin + {R max Rmin )*(x + l)/2, (53) 

(iii) MAPNO=3 (Ref. 0, 0): 

( \ P i^tan(f(l + x)) 

y[X) = Rmin H p ; — ; — , (54) 

1+ „ \ tan f(l + a: ' 1 ' 

ri.Tn.ax—ti m i n V 4 \ 



where 

p _ Rmed — Rmin ( p R n /--n 

- n 'd — p p y^max J^min)- \°° J 

Because of its flexibility, we use the tangent mapping. According to the tangent mapping, 

y(-l) = Rmin, 2/(0) = Rmed, 2/(1) = Rmax- (56) 

Therefore, one can safely control the range (R m in, Rmax) and the distribution (Rmed) of 
the grid points. With this discretization procedure, continuous integral equations are 
transformed into nonsingular matrix equations. 
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6 Description of input parameters and input distri- 
bution 



6.1 Input parameters 



NF 


Number of quark flavors 


NGRID 


Number of grid points 


PLQCD 


^qcd 


Q2I 


Initial Q 2 


Q2F 


Final Q 2 


RMIN 


Smallest grid point 


RMED 


Median of grid poiunts 


RMAX 


The maximum grid point 


LN 


The highest degree of Laguerre polynomial included 


IFLAV 


1 


u quark 




2 


d quark 


2 


s quark 


MAPNO 


1 


type-1 mapping 


2 


linear mapping 


3 


tangent mapping 


IALPH 





leading order (LO) in a s 


1 


next-to-leading order (NLO) in a s 



6.2 Arrays 



P(I) 


Grid points y = P(I), where RMIN < P(I) < RMAX 


WP(I) 


Weights 


PT(I) 


Grid points t = PT(I), where < t < T = -2/(3 \n(a s (q 2 )/a s (q%)) 


WT(I) 


Weights 


PN(LN,IE,JE) 


elements of the LN'th Laguerre moment of LO kernel [P]2x2 


RN(LN,IE,JE) 


elements of the LN'th Laguerre moment of NLO kernel [i?]2 X 2 


SPN(I,IE,JE) 


PN(I,IE,JE)- PN(I-1,IE,JE) 


E1(J,K) 


PN(0,J,K)/A 


E2(J,K) 


-PN(0,J,K)/A 


SA(K,N,IE,JE) 


[°n]2x2 


SB(K,N,IE,JE) 


[^„]2x2 


A(K,N,IE,JE) 


[^2X2 


B(K,N,IE,JE) 


K]2X2 


ENT(I,IE,JE) 


[^(0)2x2 


PHI0(Y,IE) 


Initial distributions xAY,q, xAGq respectively for IE = 1,2 


PHI0N(N,IE) 


Laguerre moments of initial distribution. 


PHIT(I,IE) 


Evolved distributions iAS, xAG respectively for IE =1,2 
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The codes provided are stored in three directories named after the polarization types; 
namely transverse, longitudinal and unpolarized. In each directory there are exacutable 
files that can be used to compile the fortran codes (executable file "comp" ) and run them 
(executable file "run"). The transverse polarization case is the simplest. For transverse 
polarization one has only the nonsinglet equation. In the longitudinally polarized and 
unpolarized cases both singlet and nonsinglet equations are solved. The solutions of the 
singlet and nonsinglet equations are provided by separate codes. For the longitudinally 
polarized case the nonsinglet equation is solved by "nslong.f " , while the singlet equation 
is solved by "snlong.f". Both codes use the same input file, which is called "INPUT". 
Nonsinglet and singlet equations work in coordination with each other. Upon execution 
of the "run" command, first the nonsinglet equation is solved, and the output is written 
into data files. Next the nonsinglet equation is solved. The data produced by the solution 
of the singlet equation is read by the code which solves the nonsinglet equation. The user 
is expected to prepare the "INPUT" file, compile the codes by using command "comp" 
and then run them using the command "run" . The procedure for the unpolarized case 
is the same, with the only difference being the names of the files which are "nsunpol.f" 
and "snunpol.f". Next, we give the detailed descriptions of the programs "sn*.f" 

7 SNLONG and SNUNPOL 

Since the only difference between these programs are the kernels and the initial condi- 
tions, the explanations provided below equally apply to both. 

7.1 Main Program 

The main program starts by calling the INPUT subroutine. This subroutine reads 12 
input parameters described in Table |6TI| . Next, the grid points and the weights are 
calculated by calling the subroutine GAUS, and they are mapped into the desired interval. 
The mapping is done by the MAP subroutine. 

Then, in order to calculate the Laguerre moments of the kernel and input distribu- 
tions, the subroutine XPL is called. The factors A^, B% which are used in the construction 
of the leading order evolution operator E® are calculated in a subroutine called XBKN. 
In the next step, the XEN subroutine is called to calculate the evolution operator E n . 
Finally, the evolution of the initial distribution is performed in the VQD subroutine. 
The results for the evolved distributions AE(x), and AG(x) are written, respectively, 
into the files "dltsig.dat", and "dltglu.dat". Later, the result for the S distribution is 
read by the program that solves the nonsinglet equation. 

7.2 Subroutine INPUT 

In this subroutine 12 input parameters are read. In addition to reading these parameters, 
various constants to be used throughout the program are also defined in the INPUT 
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subroutine. These constants are PLQCD = \ QCD , PI = tt, Z3 = ( 3 , CF, CG, CA, TR, 
BTO = fa, BT1 = Pi, T = -2/[3 ln(a s (q 2 )/a s (qi)), and XLMBD = A. 

7.3 Subroutine XPL 

In this subroutine we calculate the Laguerre moments of the initial distributions 
PHI0(Y,IE) and kernels P and R = P 1 - (3 1 /(2f3 )P . In PHI0(Y,IE), 

Y = ln(l/X) = P(J), (1) 

and the IE index is used to label the initial distributions for xAE(x)(IE=l), and 
xAG(x)(lE=2). We use variable Y = [0, inf) rather than X = [0,1] in calculating 
the Laguerre moments. The results for the Laguerre moments of the initial distribu- 
tions are stored in PHI0N(I,IE) where IE = 1,2 corresponds to the Laguerre moments 
of xAE, xAG respectively, and / refers to the order of Laguerre moments. In order to 
calculate the Laguerre moments of the kernel, we express it in the following form 

P = Pr + P \ + P S S(l - X). (2) 
I 1 _ X ) + 

The singular part P s is regulated according to the "+" regularization prescription, while 
the regular part P r is directly integrated. The delta function part requires no numerical 
integration. Therefore, the contribution of this piece is trivially added after performing 
the numerical integrals for P r and P s . The kernels are stored in the arrays P0R(Y,IE,JE), 
P0S(Y,IE,JE), P1R(Y,IE,JE), P1S(Y,IE,JE), where "R" and "S" stand for the "regular" 
and "singular" parts of the kernels, and IE = 1,2, JE = 1,2 are matrix indices. Ac- 
cording to our convention, the distributions are introduced by define 2x1 vector array 
as (xAE, xAG). In this subroutine we also define the arrays SPN(I,IE,JE), E1(J,K), and 
E2(J,K). 

7.4 Subroutine XBKN 

In this subroutine we construct the arrays SA(K,N,IE,JE), SB(K,N,IE,JE), 
A(K,N,IE,JE), B(K,N,IE,JE) using nested loops. 

7.5 Subroutine XEN 

This is the subroutine where the Laguerre moments of the evolution operator, 
ENT(N,IE,JE), are constructed. In this subroutine the functions E0N(N,T,IE,JE) and 
E1N(N,T,IE,JE) are called. These functions represent the 0th and 1st order contribu- 
tions to the Laguerre moments. 
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7.6 Function EON(N,T,IE,JE) 



Calculates EON(N,T,IE,JE) using the previously calculated A(K,N,IE,JE), and 
B(K,N,IE,JE) arrays. Notice that we have defined EON(N,T,IE,JE) as a function rather 
than a subroutine. The reason for this is the following: EON(N,XT,IE,JE) appears in- 
side an integral in the calculation of E1N(N,T,IE,JE), where XT(0 < XT < T) is the 
integration variable. Therefore, on needs to know EON for all possible values of XT. 



7.7 Function E1N(N,T,IE,JE) 

Calculates E1N. The set of grid points for the XT integral is PT(I), and the weights are 
WT(I). Grid points, which was determined in the mained program using Gauss-Legendre 
subroutines, are chosen such that < PT(I) < T. 



7.8 Subroutine VQD 

Evaluates the final result PHIT(I,IE) for evolved distributions. Here, "I" represents the 
x (Bjorken) value, where x = &• — P(I)), and IE as usual refers to two different parton 
distributions, that is PHIT(1, 1) = xG(x), and PHIT(I,2) = xE(x) (IE=2). 



7.9 Functions XL AG, FCTRL, S2, SI, ST 

XLAG(N,Y) computes the Laguerre polynomial of order N at Y. FCTRL(N) computes 
N!. The result is given in double precision. S2(X) and Sl(x) are respectively the functions 
Si(x), and S 2 {x). ST(x) represents S(x). 



7.10 Other functions 

We also define various simple functions which are used throughout the program. 
ALPHAS(Q 2 ) is a s (Q 2 ). GS(X,A,B,C,D,E,F) is the Gehrmann-Stirling ansatz "A" [|0| 
for the polarized parton distributions. In addition, for convenience in typing in the ker- 
nel for longitudinal parton distributions, we define PF(X), PGS(X), PGR(X), PNFS(X), 
PNFR(X), PA(X), FQQ(X), F1QG(X), F2QG(X), F1GQ(X), F2GQ(X), F3GQ(X), 
FIGGS(X), FIGGR(X), F2GG(X), F3GGS(X), F3GGR(X). 



7.11 subroutines GAUSS, LEGENDRE and MAP 

The LEGEND(X,L,PSUBL) subroutine computes Legendre polynomials of argument x 
from to order L. The GAUS(Y,WY,N) subroutine determines N gaussian points(in 
the vector Y) and N weights(in the vector WY). For this routine N must be even 
and not greater than 100. Since the points and weights are symmetric about 
zero, only half are stored. The original version of this subroutine was written 



by S. Cotanch [19], at the Univ. of Pittsburgh in the years 1974-75. The 
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MAP(Y,WY,N,MTYPE,RMIN,RMAX,RMED) subroutine maps the initial set of N grid 
points Y(I) and weights WY(I) to the desired interval (RMIN,RMAX). As explained in 
detail earlier, the mapping type MTYPE=3 allows one to control the median(RMED) 
of the distribution. For the Y variable that we use in calculating the Laguerre moments, 
we have used RMIN=0.0D0, RMAX=1.0D2, and RMED=5.0dO Since the Laguerre mo- 
ment integrals are damped by exponential factors e" y , RMAX=1.0D2 was a large enough 
cutoff for the integral. 

8 NSLONG and NSUNPOL 

The structure of the codes NSLONG and NSUNPOL are the same as that of SNLONG 
and SNUNPOL. All subroutine names and their functions are identicle. The only major 
difference comes from the fact that in the nonsinglet case one no longer has a matrix 
equation. Rather, there are two uncoupled equations. The nonsinglet codes read the 
xY<(x) results which are produced by the singlet codes. 

9 Running the code 

In order to get reliable results for 0.001 < x < 1 we have used approximately 30 Laguerre 
polynomials. The number of grid points used was N = 60. The nonsinglet codes run in 
about 1 minute. The singlet and nonsinglet codes together run in about 20 minutes. At 
very small x (x < 0.001) the Laguerre polynomials diverge, and therefore they are not 
a convenient basis to use in the very small x region. However, for 0.001 < x < 1, the 
results are stable for a reasonable number of Laguerre polynomials. 

10 Conclusions 

We have shown that the Laguerre expansion is a significant tool in the analysis of the 
QCD evolution equations from the numerical side. These advantages include both short 
running times in the actual implementation of the evolution and the possibility to have 
well defined recursion relations. The Laguerre algorithm is a very powerful way to address 
efficiently these problems. 

We remark that polynomial expansions are going to be of wide use in the analysis of 
more general parton distributions -such as the non-forward or the double distributions- 
which have been introduced in the recent literature on Compton processes. Here we have 
just began our tour on the analysis of QCD renormalization group equations and their 
solutions by finite step integrations. We hope to return to the study of the extension of 
these methods to the non-diagonal partonic evolution in the near future. 
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Figure 1: Evolution of the unpolarized valence quark distributions iu v versus x for 
various Q 2 values 

11 Appendix. List of all the NLO polarized kernels 

The references for the unpolarized NLO kernels are & |3j. The longitudinally polarized 
kernels can be found either as anomalous dimensions or splitting functions in [JT3j, [14] 
The Perturbative expansion of the kernels is 

AP«(*,a„) = (g) A4 0) (x) + (g) Aif (*) + ... (3) 

where ij are flavour indices. The non-singlet and the singlet polarized LO kernels, 
respectively, are given by 

a4°!(*) = P^ix) = C F -l-x + - x)\ 
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Figure 2: Unpolarized xd v versus x for various Q 2 values 



APW(x) = APJS(x 



AP£\x) = 2n f T R (2x-l) 
AP®(x) = C F {2-x) 

The "+" distributions are defined by 



2x + l] 



The non-singlet NLO kernels are given by 



o 1 — x 



Ap(l) _ ^2 



P F (x) =F Pa(x) + 5(1 - x) ( | - y + 6C(3) 
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Figure 3: Unpolarized xs distribution versus x for various Q 2 values 
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Figure 4: Unpolarized gluon distribution xG versus x for various Q 2 values 



where 

1 7T 2 

S 2 (x) = -2Li 2 (-x) - 2 In x ln(l + x) + - In 2 x . 

2 6 

The NLO polarized singlet kernels are given by 

AP$(x) = AP« + + 2C F T R n f AF gq 

AP$(x) = C F n f T R AF${x) + C G n f T R AF$ {x) 

AP^(x) = C F n f T R AF^(x) + C F F$(x) + C F C G AF^{x) 

AP^(x) = -C G T R n f AFM(x) - C F T R n f AF$(x) + C G AF$(x) 

where C G = Ca = N c and 
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Figure 5: Unpolarized xE = x(u + d + s) versus x for various Q 2 values 



AF gg (x) = 1 — x — (1 — 3x) lnx — (1 + x) In 2 x 

-22 + 27x - 9 lnx + 8(1 - x) ln(l - x) 



AF«(x 



x 



+5p qg {x) 



2 In 2 x(l — x) — 4 ln(l — x) lnx + In 2 x n z 



2(12 - liar) - 8(1 - x) ln(l - x) + 2(1 + 8x) lnx 



ln 2 (l -x) 



Tl- 



5p qg (x) — 25*2(3;) — 31n 2 x 5p, 



qg v 



-~(x + 4) - ln(l-x) 

-- - ~(4 - x) lnx - 5p gq (-x) ln(l - x) + 



-4 - ln 2 (l - x) + -ln 2 x 



Sp qg {x) 



20 



0.4 



Q 2 =4 (GeV) 2 




Figure 6: Longitudinally polarized xAu v , plotted versus x, for various Q 2 values 
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Figure 7: Longitudinally polarized a;Ad v versus x for various Q 2 values 



Sp qg (x) = 2x - 1 

SPgg(x) = \ 2X + 1. 
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The unpolarized kernels, to LO are given by 



,(0) 



qq,NS — ^ NS 



AP, 
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for the non-singlet sector, and by 



Pqq\ x ) — Pqq,NS 

P$(x) = 2T R n f (x 2 + (l-x) : 
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Figure 8: Longitudinally polarized xAs versus x for various Q 2 values 



l + (l-x) 



X 



in the singlet sector. 

The NLO unpolarized non-singlet and singlet kernels are given by 



P^(x)=AP^(x) 



>(1)T, 



and 



pW(x) = P$l + + 2C F T R n f F qq 

P«(x) = C F n f T R F$(x) + C G n f T R F$( 
3(1 

99 
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P£>{x) = C F n f T R F£\x) + C 2 F F$(x) + C F C G F$(x) 
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Figure 9: Longitudinally polarized xAG versus x, shown for various Q 2 values 



P«(x) = C F T R n f F$(x) + C G T R n f F${x) + +(%F®(x) 



(15) 



20 56 8 

F„„(x) = 2 + 6x x 2 + (1 + 5x H — re 2 ) lnx — (1 + x) In 2 x 

9x 9 3 

= 4 - 9x - (1 - Ax) lnx - (1 - 2x) In 2 x + 41n(l - x) 

2 



+ 



, o / 1 — X \ , , / 1 — X 

2 In 2 ( ) -41n 

x 



-7T Z + 10 



, 2) 182 14 40 .136 38., A . M , , . . 2 

F„ (2) = 1 x H h ( x ) In x - 4 ln(l - x) - (2 + 8x) In 2 x 

9 9 9x v 3 3 7 v ; v j 



19 



+ 



44 



In 2 x + y In x - 2 ln 2 (l - x) + 4 ln(l - x) + — 

+2p <?9 (-x)S , 2 (x) 



7T 2 218 



9 



Pqg(x) 



24 



X 



0.20 



0.10 



0.00 



-0.10 
0.01 



— Q 2 =4 (GeV) 2 
--- Q 2 =50 (GeV) 2 
Q 2 =200 (GeV) 2 




0.21 



0.41 



0.61 



0.81 



Figure 10: Longitudinally polarized xAE versus x, shown for various Q 2 values 
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Figure 11: Transversely polarized xAtu v versus x, shown for various Q 2 values 
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Figure 12: Transversely polarized xA-rd v versus x, shown for various Q 2 values 



The kernel for the transverse polarization is written as 
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The NLO corrections are given by [IT 
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Figure 13: Transversely polarized xAts, versus x, shown for various Q 2 values 
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